###This file is written by Gunnar Waterstraat###
###Clinic of neurology, Charite CBF###

#################Zero-phase-shift filter from Scipy-Cookbook########################
####################################################################################
from numpy import vstack, hstack, eye, ones, zeros, linalg, \
newaxis, r_, flipud, convolve, matrix, array
from scipy.signal import lfilter

def _lfilter_zi(b,a):
    #compute the zi state from the filter parameters. see [Gust96].
    #Based on:
    # [Gust96] Fredrik Gustafsson, Determining the initial states in forward-backward 
    # filtering, IEEE Transactions on Signal Processing, pp. 988--992, April 1996, 
    # Volume 44, Issue 4
    n=max(len(a),len(b))
    zin = (eye(n-1) - hstack( (-a[1:n,newaxis],vstack((eye(n-2), zeros(n-2))))))
    zid=  b[1:n] - a[1:n]*b[0]
    zi_matrix=linalg.inv(zin)*(matrix(zid).transpose())
    zi_return=[]
    #convert the result into a regular array (not a matrix)
    for i in range(len(zi_matrix)):
      zi_return.append(float(zi_matrix[i][0]))
    return array(zi_return)

def _filtfilt(b,a,x):
    #For now only accepting 1d arrays
    ntaps=max(len(a),len(b))
    edge=ntaps*3
    if x.ndim != 1:
        raise ValueError, "Filiflit is only accepting 1 dimension arrays."
    #x must be bigger than edge
    if x.size < edge:
        raise ValueError, "Input vector needs to be bigger than 3 * max(len(a),len(b)."
    if len(a) < ntaps:
        a=r_[a,zeros(len(b)-len(a))]
    if len(b) < ntaps:
        b=r_[b,zeros(len(a)-len(b))]
    zi=_lfilter_zi(b,a)
    #Grow the signal to have edges for stabilizing 
    #the filter with inverted replicas of the signal
    s=r_[2*x[0]-x[edge:1:-1],x,2*x[-1]-x[-1:-edge:-1]]
    #in the case of one go we only need one of the extrems 
    # both are needed for filtfilt
    (y,zf)=lfilter(b,a,s,-1,zi*s[0])
    (y,zf)=lfilter(b,a,flipud(y),-1,zi*y[-1])
    return flipud(y[edge-1:-edge+1])
####################################################################################
####################################################################################

def GetMarkerPosFromData(marker, smooth=2):
    """Gets the marker position from the trigger channel
    GetMarkerPosFromData(marker)
  
    input:
    -- marker - one-dimensional array with trigger channel
                    each crossing of the mean value of that channel
                    is interpreted as a trigger
    -- smooth - integer: smoothing data by calculating sliding average of 'smooth'
    following points
    output:
    -- marker - one-dimensional array containing trigger positions
    """
    from numpy import vstack, abs
    from scipy.signal import detrend
    from scipy import where
    marker = detrend(marker.reshape([1,-1]))
    marker_temp = marker[0,0:marker.shape[1]-smooth]
    for i in xrange(1,smooth,1):
        marker_temp = vstack([marker_temp,marker[0,i:marker.shape[1]-smooth+i]])
    marker = marker_temp.mean(axis=0)
    cross_mult = marker[:-1] * marker[1:]
    PosToNeg = where((cross_mult < 0) & (marker[:-1] > 0),-1,0)
    NegToPos = where((cross_mult < 0) & (marker[:-1] < 0),1,0)
    marker = PosToNeg + NegToPos
    return where(abs(marker)==1)[0]

def _GetBPFilterCoeffs(data,lp=400,hp=1200,sf=5000,width=30,gpass=5,gstop=10):
    """Helper function for BPfilter
    """
    from scipy.signal import buttord, butter
    ny = sf/2.
    lp /= ny
    hp /= ny
    width /= ny
    [order,wn] = buttord(wp=[lp,hp],ws=[lp-width,hp+width],gpass=gpass,gstop=gstop,analog=0)
    buttcoeffs = butter(order,wn, btype='bandpass')
    return buttcoeffs[0],buttcoeffs[1]

def BPFilter(data,lp=400,hp=1200,sf=5000,width=30,gpass=5,gstop=10):
    """Bandpass filter the data,
    BPFilter(data,lp=400,hp=1200,sf=5000,width=30,gpass=5,gstop=10)

    input:
    -- data - one-dimensional array containing the channel to be bp-filtered
    -- lp   - Lower border of the frequency band
    -- hp   - Higher border of the frequency band
    -- sf   - sampling frequency
    -- width - width of the filter slope
    -- gpass - maximal attenuation of the signal in the passband in dB
    -- gstop - minimal attenuation of the siganl in the stopband in dB
    output:
    -- filt - filtered array
    """
    coeff1,coeff2 = _GetBPFilterCoeffs(data,lp=lp,hp=hp,sf=sf,width=width,gpass=gpass,gstop=gstop)
    return _filtfilt(coeff1,coeff2,data)

def CalculateRMS(data):
    """
    Calculate rms value of the input data along axis 0
    Calculaterms(data)

    input:
    -- data - input data
    output:
    -- rms value along axis 0
    """
    from scipy import mean, sqrt, shape, ndim
    data = data-mean(data,axis=0)
    rms=sqrt((data**2).sum(axis=0)/shape(data)[0])
    if ndim(data)==3:
        rms=rms.swapaxes(0,1)
    return rms

###Code for Calculating spherical spline interpolation and CSD maps########
###Algorithm from  Perrin et al., Electroenceph Clin Neurophysiol##########
###1989,72:184-187, Corrigenda in 1990, 76: 565-566########################
###########################################################################

def MeshCircle(d_samples=100):
    """
    Calculate a meshgrid containing the points of a circle (center (0,0), radius 1)
    MeshCircle(d_samples=100)

    input:
    -- d_samples - number of points along the diameter
    output
    -- coords - 2D array containing the coordinates in rows
              - 1st column: x
              - 2nd column: y
    """
    from numpy import array, linspace, sqrt, ones, vstack, hstack
    from scipy import where
    coords = array([0,1])
    y = linspace(-1.,1.,d_samples,endpoint=True)
    for el in y[1:-1]:
        endpoint = sqrt(1.**2 - el**2)
        coords = vstack([coords,array([-1*endpoint,el])])
        temp_x = y[where((y > (-1 * endpoint)) & (y < endpoint))[0]]
        coords =\
        vstack([coords,vstack([temp_x,ones(len(temp_x))*el]).swapaxes(0,1)])
        coords = vstack([coords,array([endpoint,el])])
    coords = vstack([coords,array([0,-1.])])
    return coords

def MeshCircleOnSphere(coords):
    """
    Calculate a 3D meshgrid on a sphere from a 2D meshgrid on a circle
    MeshCircleOnSphere(coords)

    input:
    -- coords - 2D array containing the coordinates in rows
              - 1st column: x
              - 2nd column: y
    output
    -- sphere_coords - 2D array containing the coordinates in rows
                     - 1st column: x
                     - 2nd column: y
                     - 3rd column: z
    """
    from numpy import sqrt, sin, empty
    sphere_coords = empty([coords.shape[0],3],dtype=float)
    theta = sqrt(coords[:,0]**2 + coords[:,1]**2)
    k = sin(theta) / theta
    sphere_coords[:,0] = k*coords[:,0]
    sphere_coords[:,1] = k*coords[:,1]
    sphere_coords[:,2] = sqrt(1 - sphere_coords[:,0]**2 - sphere_coords[:,1]**2)
    return sphere_coords
      
def GetStandardCoordinates(elec_names,fname):
    """
    Read Electrode Coordinates from tab-seperated text file (standard: sphere_1005.txt)
    GetStandardCoordinates(elec_names)

    input:
    -- elec_names - List of Electrode names
    output:
    -- coords - 2D array containing the coordinates in rows
                    - 1st column: x
                    - 2nd column: y
                    - 3rd column: z
    """
    from csv import reader
    from numpy import array, abs, nan, matrix
    filereader = reader(open(fname), delimiter='\t')
    coords = dict([(row[0],row[1:4]) for row in filereader])
    RP = coords['RP']
    NA = coords['Na']
    CZ = coords['Cz']
    TransMat = matrix(array([RP, NA, CZ], dtype=float).T)
    coords_result = []
    for item in elec_names:
        if item in coords:
              temp = list(array(TransMat.I * matrix(array(coords[item], \
              dtype=float)).T).flatten())
              coords_result.append(temp)
        else: coords_result.append(array([nan,nan,nan]))
    coords = array(coords_result,dtype=float)
    #if fname == 'sphere_1005.txt':
    #      coords[:,:2] *=10
    #      coords[:,2] = abs(10 * coords[:,2])
    return coords

def GetChannelNames(fname):
    """
    Read Names of Electrodes from text file
    GetChannelNames(fname)

    input:
    -- fname - filename of textfile containing electrode-names
    output:
    -- elecnames -  List of electrode names
    """
    from csv import reader
    from numpy import array
    filereader = reader(open(fname), delimiter='\t', quotechar='"')
    elecnames = dict([(row[0],int(row[1])) for row in filereader])
    from operator import itemgetter
    elecnames = sorted(elecnames.iteritems(), key=itemgetter(1))
    return [item[0] for item in elecnames]

def _CalcCosBetwElec(elec1,elec2):
    """
    Calculate Cosine between to Electrodes
    _CalcCosBetwElec(elec1,elec2)

    input:
    --elec1 - Cartesion Coordinates of Electrode1
                    - 1st column: x
                    - 2nd column: y
                    - 3rd column: z
    --elec2 - Cartesion Coordinates of Electrode2
                    - 1st column: x
                    - 2nd column: y
                    - 3rd column: z
    output:
    -- cos: Cosine between both electrodes
    """
    from numpy import isnan
    cos = 1 - ((elec1-elec2)**2).sum() / 2.
    if isnan(cos): 1/0
    return cos

def _CalcP(cos, n=7):
    """
    Calculate the first n terms of Legendre Polynomial
    _CalcP(cos, n)

    input:
    -- cos - Cosine between two electrodes
    -- n - Number of terms to be Calculated
    output:
    -- p - array containing Legendre polynomial from term 0 at p[0] to n at p[n]
    """
    from numpy import array, hstack
    p = array([1,cos]).astype(float)
    for i in xrange(1,n,1):
        p = hstack([p, ((2*i + 1) * cos * p[i] - i * p[i-1])/ float(i+1)])
    return p

def _CalcG(cos, n=7, m=4):
    from numpy import pi, arange
    p = _CalcP(cos,n)
    n = arange(1,n+1,1)
    g = 1/float(4*pi) * ((2*n + 1)/(n**m *(n+1)**m).astype(float) * p[1:]).sum()
    return g

def _CalcH(cos, n=20, m=4):
    from numpy import pi, arange
    p = _CalcP(cos,n)
    n = arange(1,n+1,1)
    h = (1/float(4*pi)) * ((2*n + 1)/(n**(m-1) * (n+1)**(m-1)).astype(float) * p[1:]).sum()
    h1 = (1 / float(4*pi)) * (((2*n + 1) / (n**(m-1) * \
    (n+1)**(m-1)).astype(float)) * p[1:]).sum()
    return h

def _CalcGij(InterpCoords,RealCoords, n=7, m=4, addzeros=False):
    from numpy import array, empty
    Gij = empty([RealCoords.shape[0]+1, InterpCoords.shape[0]+1])
    for i in xrange(InterpCoords.shape[0]):
        Gij[:-1,i] = \
        array([_CalcG(_CalcCosBetwElec(InterpCoords[i],RealCoords[j]), n=n, m=m) \
        for j in xrange(RealCoords.shape[0])])
    if addzeros == True:
        Gij[:,-1] = hstack([ones([RealCoords.shape[0]]),0])
        Gij[-1,:] = hstack([ones([InterpCoords.shape[0]]),0])
        return Gij
    else:
        return Gij[:-1,:-1]

def _CalcHij(InterpCoords,RealCoords, n=20, m=4):
    from numpy import array, empty
    Hij = empty([RealCoords.shape[0], InterpCoords.shape[0]])
    for i in xrange(InterpCoords.shape[0]):
        Hij[:,i] = \
        array([_CalcH(_CalcCosBetwElec(InterpCoords[i],RealCoords[j]), n=n, m=m) \
        for j in xrange(RealCoords.shape[0])])
    return Hij

def CalcHemSplineMap(RealCoords, DP, diameter_samples=30, type='2d', n=7, m=4):
    """Interpolate Points on a sphere using spherical splines
    CalcHemSplineMap(RealCoords, DP, diameter_samples, type, n, m)

    input:
    -- RealCoords - (x,y,z) cartesian coordinates of physical electrodes
    -- DP - array containing values a electrodes in same order as in 'RealCoords'
    -- diameter_samples - number of points along the diameter of the scalp to
    interpolate
    -- type - Dimension of ouput electrode postions, yet only '2d' implemented
    -- n - How many terms of the Legendre-Polynomial should be calculated, defaults
    to 7
    -- m - order of the spherical spline interpolation, defaults to 4
    output:
    -- InterpCoords - array containing (x,y,z) coordinates of interpolated points
    -- pot - array containing values at interpolated points
    """
    from numpy import hstack, array, dot, sum
    from scipy.linalg import pinv2
    InterpCoords_2d = MeshCircle(d_samples=diameter_samples)
    InterpCoords_3d = MeshCircleOnSphere(InterpCoords_2d)
    Gij =  _CalcGij(RealCoords, RealCoords, n=20, m=m, addzeros=True)
    DP = hstack([DP,array([0])])
    DP = DP.reshape([-1,1])
    c = dot(pinv2(Gij), DP)
    g =  _CalcGij(InterpCoords_3d, RealCoords, n=n, m=m)
    pot = c[-1,0] + sum(c[:-1,:]*g, axis=0)
    return InterpCoords_2d, pot

def CalcCSDSplineMap(RealCoords, DP, diameter_samples=30, type='2d', n=20, m=4):
    """Calculate current source densities on a sphere using spherical splines
    CalcCSDSplineMap(RealCoords, DP, diameter_samples, type, n, m)

    input:
    -- RealCoords - (x,y,z) cartesian coordinates of physical electrodes
    -- DP - array containing values a electrodes in same order as in 'RealCoords'
    -- diameter_samples - number of points along the diameter of the scalp to
    interpolate
    -- type - Dimension of ouput electrode postions, yet only '2d' implemented
    -- n - How many terms of the Legendre-Polynomial should be calculated, defaults
    to 7
    -- m - order of the spherical spline interpolation, defaults to 4
    output:
    -- InterpCoords - array containing (x,y,z) coordinates of interpolated points
    -- pot - array containing CSD-values at interpolated points
    """
    from numpy import hstack,array, dot, sum
    from scipy.linalg import pinv2
    InterpCoords_2d = MeshCircle(d_samples=diameter_samples)
    InterpCoords_3d = MeshCircleOnSphere(InterpCoords_2d)
    Gij =  _CalcGij(RealCoords, RealCoords, n=20, m=m, addzeros=True)
    DP = hstack([DP,array([0])])
    DP = DP.reshape([-1,1])
    c = dot(pinv2(Gij), DP)
    h =  _CalcHij(InterpCoords_3d, RealCoords, n=n, m=m)
    pot = sum(c[:-1,:]*h, axis=0)
    return InterpCoords_2d, pot

def MeshScalpMap(coords,pot,gridsize=200):
    """Helper function to plot irregularly spaced data with pylab.pcolor
    MeshScalpMap(coords,pot,gridsize)

    imput:
    -- coords - array containing 2D-(x,y)-coordinates of points
    -- pot - value a points in same order as 'coords'
    -- gridsize - number of points to be plotted, output will have size
    gridsizexgridsize
    output:
    -- xi, yi, zi - x coordinates, y coordinates, values a grid position
    """
    from matplotlib.mlab import griddata
    from numpy import linspace
    xi = linspace(-1,1,gridsize, endpoint=True)
    yi = linspace(-1,1,gridsize, endpoint=True)
    zi = griddata(coords[:,0],coords[:,1],pot,xi,yi)
    return xi, yi, zi

######################################################################################
######################################################################################

####################Perform Correlation Analysis######################################
######################################################################################
def CalcSpearmanAutocorr(data, n=10):
    """
    Calculate Spearman Autocorrelation Coefficient
    Usage: CalcSpearmanAutocorr(data, n=10)

    --input:
        data: array of input values
        n: amount of shift up to which autocorrelation should be calculated
    --output:
        array with shape (2,n[,n_channels])
        -first axis:
            -first element is spearman r
            -second element is two-tailed p-value
    """
    from numpy import array
    import scipy.stats as Sst
    if data.ndim == 1:
        result = array([Sst.spearmanr(data[0:-n],data[k:k-n]) for k in \
             xrange(n)]).swapaxes(0,1)
    else:
        result = []
        for ch in xrange(data.shape[1]):
            result.append(array([Sst.spearmanr(data[0:-n,ch],data[k:k-n,ch]) for k in \
                 xrange(n)]))
        result = array(result).swapaxes(0,2)
    return result

def CalcSpearmanCorr(data1, data2, n=10):
    """
    Calculate Spearman Correlation Coefficient
    Usage: CalcSpearmanCorr(data1, data2, n=10)

    --input:
        data1: array of first input values
        data2: array of second input values
        n: amount of shift up to which correlation should be calculated
    --output:
        array with shape (2,n[,n_channels])
        -first axis:
            -first element is spearman r
            -second element is two-tailed p-value
    """
    from numpy import array
    import scipy.stats as Sst
    if data1.ndim == 1:
        result = array([Sst.spearmanr(data1[0:-n],data2[k:k-n]) for k in \
             xrange(n)]).swapaxes(0,1)
    else:
        result = []
        for ch in xrange(data1.shape[1]):
            result.append(array([Sst.spearmanr(data1[0:-n,ch],data2[k:k-n,ch]) for k in \
                 xrange(n)]))
        result = array(result).swapaxes(0,2)
    return result

def CalcSpearmanCorrShuffledSides(data1, data2, n=10):
    """
    Calculate Spearman Correlation Coefficient with elements of 'data1' and
    'data2' shuffled randomly from side to side
    Usage: CalcSpearmanCorr(data1, data2, n=10)

    --input:
        data1: array of first input values
        data2: array of second input values
        n: amount of shift up to which correlation should be calculated
    --output:
        array with shape (2,n[,n_channels])
        -first axis:
            -first element is spearman r
            -second element is two-tailed p-value
    """
    from numpy import array, random, newaxis
    from scipy import where
    data1 -= data1.mean(axis=0)[newaxis,:]
    data2 -= data2.mean(axis=0)[newaxis,:]
    data1 /= data1.std(axis=0)[newaxis,:]
    data2 /= data2.std(axis=0)[newaxis,:]
    indices = random.random_integers(0,1,data1.shape)
    data = where(indices==0,data1,data2)
    return CalcSpearmanAutocorr(data, n=n)

def CalcSpearmanCorrShuffledTrials(data1, data2, n=10):
    """
    Calculate Spearman Correlation Coefficient with elements of 'data1' and
    'data2' shuffled randomly trial-by-trial
    Usage: CalcSpearmanCorr(data1, data2, n=10)

    --input:
        data1: array of first input values
        data2: array of second input values
        n: amount of shift up to which correlation should be calculated
    --output:
        array with shape (2,n[,n_channels])
        -first axis:
            -first element is spearman r
            -second element is two-tailed p-value
    """
    from numpy import array, random, newaxis
    from scipy import where
    data1 -= data1.mean(axis=0)[newaxis,:]
    data2 -= data2.mean(axis=0)[newaxis,:]
    data1 /= data1.std(axis=0)[newaxis,:]
    data2 /= data2.std(axis=0)[newaxis,:]
    indices = (random.random(data1.shape).argsort(axis=0)).astype(int)
    data1 = array([data1[indices[:,k],k] for k in \
                   xrange(indices.shape[1])]).swapaxes(0,1)
    data2 = array([data2[indices[:,k],k] for k in \
                   xrange(indices.shape[1])]).swapaxes(0,1)
    return CalcSpearmanCorr(data1, data2, n=n)
######################################################################################
######################################################################################
      
def ReadBinary(fname,num_channels,channels='all',readnum_dp='all',data_type='float4', buffermem=512):
    from os.path import getsize
    from math import ceil
    from numpy import arange, empty
    if data_type == 'float4':
        bytenum = 4
        id = 'f'
    filesize = getsize(fname)
    data_num = filesize / bytenum / num_channels
    if channels =='all': channels = arange(num_channels)
    if readnum_dp =='all': readnum_dp = data_num
    from scipy.io.numpyio import fread #import module to read binary data
    fd = open(fname,'rb') #open file
    bytequot = int(buffermem * 1024**2 / bytenum) - (buffermem * 1024**2 / bytenum) % num_channels
    batchnum = int(ceil(num_channels*readnum_dp / float(bytequot)))
    num_dp_per_batch = bytequot / num_channels
    readnum_ch = len(channels)
    data = empty([readnum_ch,readnum_dp],dtype=id)
    if (num_channels*readnum_dp) % bytequot == 0:
        for i in xrange(batchnum):
              data_temp=fread(fd, bytequot,id).reshape([num_channels,-1],order='f')
              data[:,i*num_dp_per_batch:i*num_dp_per_batch + num_dp_per_batch] \
              = data_temp[channels,:]
        fd.close()
    else:
        for i in xrange(batchnum-1):
              data_temp=fread(fd, bytequot,id).reshape([num_channels,-1],order='f')
              data[:,i*num_dp_per_batch:i*num_dp_per_batch + num_dp_per_batch] \
              = data_temp[channels,:]
        data_temp=fread(fd,(num_channels*readnum_dp) % bytequot,id).reshape([num_channels,-1],order='f')
        data[:,(batchnum-1)*num_dp_per_batch:] = data_temp[channels,:]
        fd.close()
    print 'data read'
    return data

def InterpolationFromData(data, markers, win, mode='ms'):
    """interpolates the stimulus artefacts in wanted time window

    data=interpolationfromfile(data, markers, win, mode)
    input:
    -- data - two dimensional array
              1st dimension: data points
              2nd dimension: channel number from 0...n
    -- markers - marker positions arranged in 1d array
    -- win - time window to be interpolated arranged in array (time in ms or dp before stim,time in ms or dp after stim), only positive values allowed
    -- mode - 'ms' for win given in miliseconds or 'dp' for win given in datapoints
    output:
    -- data - two dimensional array
              1st dimension: data points
              2nd dimension: channel number from 0...n
    """
    from scipy import absolute, arange, reshape, array, ndim
    from scipy.interpolate import interp1d
    if mode=='ms':
	    win=(absolute(win)*fs/1000).astype(int)
    elif mode=='dp':
	    win=absolute(win)
    else:
	    print 'mode not understood'

    # Calculate datapoints to be interpolated
    interpolpts=[arange(i-win[0],i+win[1],1) for i in markers]
    interpolpts=reshape(interpolpts,(-1,))
    
    #set fixed data points for interpolation
    x=array((markers-win[0],markers+win[1]))
    x=reshape(x,(-1,))
    x.sort()        # y would be data[x,ch]
    
    if ndim(data)==2:
	    (lenodata, noofusedchan)=data.shape
	    for ch in xrange(noofusedchan):
		f=interp1d(x,data[x,ch].astype(float), kind='linear')
		data[interpolpts,ch]=(f(interpolpts)).astype(float)
		print 'channel '+str(ch)+' interpolated'
    elif ndim(data)==1:
	    f=interp1d(x,data[x].astype(float), kind='linear')
	    data=data.astype(int)
	    data[interpolpts]=(f(interpolpts)).astype(float)
	    print 'data interpolated'
    print 'interpolating done'
    return data

def GetTrialsFromData(data,stim,win,mode='ms'):
    """reshape data into single trials matrix:

     -- stim - indices of stimuli (in samples)
     -- data - eeg data
     -- win - data window  in ms or datapoints which should be arranged as array or list
     -- mode - 'ms' or 'dp' for unit of window, standard is 'ms'
    """ 

    from scipy import fix, array, empty, where, ndim
    if mode=='ms':
      nstim=len(stim)-2
      first=int(fix(win[0]/1000.0*fs))
      npts=int(fix(win[1]/1000.0*fs)-first)
    elif mode=='dp':
        nstim=len(stim)-2
        first=win[0]
        npts=win[1]-first
    else:
      print 'mode not understood'
    if data.ndim==2:
      (datalen,nchans)=data.shape
      (idx,)=where(array(stim)<datalen)
      idx=idx[:-1]
      trials=empty((npts,nchans,len(idx)))
      for i in idx:
          trials[:,:,i]=(data[stim[i]+first:stim[i]+first+npts,:])
    elif ndim(data)==1:
      datalen = len(data)
      (idx,) = where(array(stim)<datalen)
      idx=idx[:-1]
      trials = empty((npts,len(idx)))
      for i in idx:
          trials[:,i]=(data[stim[i]+first:stim[i]+first+npts])
    return trials

if __name__=='__main__':
    import numpy as N
    import pylab as P
    import scipy as S
    from scipy.signal import detrend
    fname = 'Channel_names.csv'
    elecnames = GetChannelNames(fname)
    elecnames = N.hstack([N.array(elecnames[:14]),N.array(elecnames[16:])])
    RealCoords = GetStandardCoordinates(elecnames) 
    channels = N.arange(128,160,1)
    wanted_channels = \
    N.array([0,1,2,3,4,5,6,7,8,9,10,11,12,13,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30])
    data = ReadBinary('Gunnar1-export.bin',160,channels) * 3.6328e+014
    data = detrend(data, axis=-1)
    marker = GetMarkerPosFromData(data[-1,:])
    data = data[:-1,:]
    marker_about = S.array([N.arange(m+500,m+515,1) for m in marker]).astype(int)
    stim_art_pos = N.array([S.where(N.abs(data[4,marker_about[k]]) == \
    N.abs(data[4,marker_about[k]]).max())[0] for k in xrange(len(marker))]).flatten()
    from scipy.stats import mode
    marker2 = marker + 500 + (stim_art_pos - mode(stim_art_pos)[0]).astype(int)
    data = InterpolationFromData(data.swapaxes(0,1), N.sort(N.hstack([marker,marker2])), \
    [0,20], mode='dp')
    t = N.linspace(0,90,450,endpoint=False)
    yfilt2 = \
    N.array([BPFilter(data[:,k],lp=500,hp=900,sf=5000,width=30,gpass=5,gstop=10) \
    for k in xrange(data.shape[1])]).swapaxes(0,1)
    bp2_av1 = GetTrialsFromData(yfilt2,marker,[0,450],mode='dp').mean(axis=2)
    bp2_av2 = GetTrialsFromData(yfilt2,marker2,[0,450],mode='dp').mean(axis=2)
